Activated sampling in complex materials at finite temperature: the 
properly-obeying-probability activation-relaxation technique 
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While the dynamics of many complex systems is dominated by activated events, there are very 
few simulation methods that take advantage of this fact. Most of these procedures are restricted 
to relatively simple systems or, as with the activation-rela:xation technique (ART), sample the 
conformation space efficiently at the cost of a correct thermodynamical description. We present 
here an extension of ART, the properly-obeying-probability ART (POP-ART), that obeys detailed 
balance and samples correctly the thermodynamic ensemble. Testing POP-ART on two model 
systems, a vacancy and an interstitial in crystalline silicon, we show that this method recovers the 
proper thermodynamical weights associated with the various accessible states and is significantly 
faster than MD in the diffusion of a vacancy below 700 K. 

PACS numbers: 5.10.-a, dynamics 5.70.-a , 66.30.-h, 68.35. Fx) 82.20.Wt simulation 
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I. INTRODUCTION 

Throughout physics, chemistry and biology, a large 
proportion of atomistic processes take place on time 
scales many orders of magnitude longer than the typi- 
cal period of atomic vibrations. These processes are out 
of reach of straightforward molecular dynamics (MD), 
which cannot exceed simulation times equivalent to the 
microseconds at best. It is not surprizing, therefore, that 
considerable effort has been devoted in the last few years, 
with some degree of success, to develop algorithms that 
overcome this hmitation 0, H II II II SM i, 113 ■ 

These methods can be separated into two classes. Ac- 
tivated methods, such as the activation-relaxation tech- 
nique (ART) P,|3 and related approaches 0,0, 111 sample 
the energy landscape of complex systems by identifying 
minima connected by minimum energy pathways. This 
family of methods is very efficient for sampling conforma- 
tions. Recently, ART was found to be the most efficient 
method for high-dimensional problems However, 
because of an unknown bias in the selection of events 
for these methods, it is not possible to ensure a proper 
thermodynamic sampling. While this is not a major lim- 
itation for sampling states or even identifying pathways, 
such as protein folding trajectories, for example, it is suf- 
ficiently severe to prevent the use of these methods to 
sample equilibrium or dynamical quantities. 

The second class of methods is based on molecular 
dynamics and corresponds to methods such as hyper- 
MD 0, temperature-assisted dynamics (TAD) Q, self- 
guided dynamics 0| and biased thermodynamics [lOj . 
Until now, the application of these methods has been 
mostly restricted to simple systems, with a limited num- 
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ber of relatively well-characterized barriers or with a well- 
defined reaction coordinate. Recently, Choudhary and 
Clancy have proposed modifications of hyper-MD that 
could allow its application to disordered materials [T^ . 
It appears, however, that applying a significant boost in 
hyper-MD could break the thermodynamical character of 
the algorithm, placing this method in the first category 
of activated methods. 



In this paper, we present an algorithm that offers cor- 
rect thermodynamical sampling without suffering from 
the usual exponential slowing down with decreasing tem- 
perature. The properly-obeying-probability ART (POP- 
ART) samples the thermodynamically relevant parts of 
phase space, hopping efficiently over high barriers sep- 
arating low-energy basins. We apply this algorithm to 
two test cases, the diffusion of a sclf-intcrstitial andj:he 
self-diffusion of a vacancy in Stillinger- Weber silicon jl^ , 
to verify the correctness of the method, and to assess its 
efficiency. 



The paper is organized as follows. We start with a 
brief discussion of limitations of standard activated meth- 
ods, such as ART. We then introduce POP-ART and 
show how it can overcome these limitations and ensure 
a proper thermodynamical sampling. The justification 
for the various steps needed to construct activated path- 
ways with detailed balance is then discussed in detail and 
a summary of the algorithm is given. The algorithm is 
tested in a study of the diffusion of an interstitial and 
a vacancy in Stillinger- Weber silicon. In Appendix A, 
we discuss a physical interpretation of the Jacobian used 
in POP-ART; in Appendix B, we present an analytical 
calculation for a simple model potential that provides 
further insights into the method; 
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II. SAMPLING THE ENERGY LANDSCAPE 
USING ACTIVATED METHODS 

The energy landscape of a system of M atoms can 
be represented as an = 3M-dimensional hypersurface, 
with the "height" indicating the value of the potential of 
the configuration at a given set of atomic coordinates. In 
a dynamical regime dominated by rare events, a system 
spends most of its time oscillating near a local energy 
minimum, hopping over an energy barrier only when a 
thermal fluctuation transfers sufficient energy onto the 
corresponding mode. Since the probability of energy 
transfer decreases exponentially with its size, the acti- 
vated trajectory will tend to cross near the lowest-energy 
point on the ridge, corresponding to a first-order saddle 
point. 

It is possible to reconstruct these trajectories, as a se- 
quence of local minima separated by transition points, 
using the activation-relaxation technique or related 
methods 0, S ■ In its latest form, called ART nou- 
veau this method works in three steps: (1) Starting 
from a local minimum, a deformation is applied until the 
lowest curvature of the Hessian matrix, given by 

OXiOXj 

becomes negative. (2) The configuration is then pushed 
along the corresponding direction while the energy in the 
perpendicular hyperplane is minimized until it reaches a 
first-order saddle point. (3) The configuration is then 
pushed slightly further, away from this saddle point, and 
its energy is minimized using a standard minimization 
technique. 

ART and similar methods have been applied with suc- 
cess to study the topology of the energy landscape and ac- 
tivated mechanisms in a wide range of materials inclu ding 
amorphous and crystalline semiconductor s 113. Iisl 0|, 
glassy materials u], atomic clusters |^ Il8j. and pro- 
teins 113, 20, 21j, (22c ,23] . A recent study has shown that 
ART compares favorably with other related techniques in 
terms of efficiency and completeness of finding activated 
mechanisms 11]. 

It is tempting to use ART approaches to study dy- 
namical trajectories, generating an on-the-fly catalog of 
events that can be selected using kinetic Monte Carlo. 
This approach was used by Henkelman and Jonsson in 
the simulation of growth on an Al (100) surface [2^ . 
and by Middleton and Wales for a binary Lennard 
Jones glass. Comparing with molecular dynamics, Mid- 
dleton and Wales observed that the kMC results with an 
on-the-fly catalog are qualitatively incorrect. It is likely 
that this discrepancy is at least partly due to an uncon- 
trolled selection bias for these types of methods, making 
it impossible to ensure a proper sampling of the barriers 
generating the catalog of kMC events. 

These limitations are fundamental and cannot be lifted 
by simply doubling or tripling the sampling. It is essen- 
tial to incorporate thermodynamics at the core of the 



algorithm. One method for achieving this is discussed in 
the next section. 



III. THE POP-ART APPROACH 

We start by separating the configuration space into two 
regions: the basin and the saddle regions. Basins are de- 
fined in a way that ensures that they contain most of the 
thermodynamical weight at a given temperature, while 
the saddle regions are visited only rarely and in passing. 
The dividing (hyper-)surface between these two regions 
is chosen to be the surface where the lowest curvature of 
the potential energy surface equals a threshold value Aq . 
The basins represent the parts of the configuration space 
where the lowest curvature has a value above Aq; they 
form a series of disconnected regions surrounding local 
minima. The saddle region is on the other side of the 
threshold and includes most other local extrema, such as 
first- and higher-order saddle points (see Fig.^. The use 
of this criterion for separating the configuration space is 
convenient as the status of any point in the configuration 
space can be decided locally, without relaxing to a nearby 
stationary point. For a given threshold, it is always pos- 
sible that the negative eigenvalue associated with a par- 
ticular saddle is higher (lower by absolute value) than the 
chosen threshold, such that it belongs to a basin. As will 
be seen from the treatment of the method below, this 
does not invalidate the algorithm but may even be used 
to one's advantage. 

Having separated the energy landscape, we define mo- 
tion in each of these regions. All the motion within the 
basin is performed with conventional MD at the desired 
temperature. Once the configuration hits the dividing 
surface, the MD is halted, the configuration is brought 
through the saddle region to a new basin at the same en- 
ergy, according to the activation rules described below, 
and the MD is resumed at the same temperature. All 
steps respect detailed balance and the overall trajectory 
samples the basins according to the proper thermody- 
namical ensemble. 

The activated part of the algorithm is composed of 
two steps: (1) the activation trajectory is first generated, 
from one basin to the other, and then (2) the free energy 
difference between the beginning and the end of this tra- 
jectory is calculated. The latter information is used for 
the accept-reject step. 

In the next subsections, we discuss these two steps in 
detail before presenting the algorithm as it is currently 
implemented. 



A. The activated trajectory 

As in ART, an activation trajectory is created by mov- 
ing along the eigenvector corresponding to the lowest 
eigenvalue of the Hessian. Unlike ART, however, there is 
no relaxation in the perpendicular hyperplane. Instead, 
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FIG. 1: Sketch of a two-dimensional energy landscape. The 
black dots denote the locations of local-energy minima. These 
minima are part of basins, bounded by a line of constant 
lowest-curvature (solid line); the percolating region surround- 
ing the basins is called the saddle region. Basin-to-basin tra- 
jectories as generated by THWART are indicated by dashed 
lines. Constrained to ensure detailed balance, the trajectories 
come back to were they started if they fail to find a boundary. 

all atoms are moved in such a way as to keep the to- 
tal potential energy constant. This is easily achieved as 
the configuration is thermalized, with roughly ksT /2 of 
available potential energy per degree of freedom above 
that in the local energy minimum. More specifically, the 
activated trajectory is generated by iterating the follow- 
ing equation: 

X,+i =X^ + ^{h^ + -I- ^ + Fi+i^ , (2) 

where, hi is the normalized eigenvector at Xi correspond- 
ing to the lowest Hessian eigenvalue, Fi is the total iV- 
dimensional force at afi, Ar is a constant that determines 
the size of the increment, and c is a multiplicative con- 
stant, chosen to project the trajectory onto the hyper- 
surface of constant potential energy. 

The orientation of is chosen initially so as to point 
in the direction of more negative curvature, i.e., away 
from the initial basin; it is updated at each step by re- 
quiring that the inner product of the local eigenvector hi 
with that at the previous step, be always positive. 

Values of h and F at point x^+i are obtained iteratively, 
i.e., initially /i^+i = hi and F^+i = Fi are used in Eq. ^ 
to obtain a value of x^+i, then values of h and F are 
calculated at the new point and inserted into Eq. |(5J to 
get the next iteration, etc. 

Unlike in ART, there is no separate relaxation stage, 
and Eq. Q is iterated across the saddle region until 
the new basin is reached, i.e., until the lowest eigen- 
value passes the threshold (from below, this time). At 
this point, the activation-relaxation phase is stopped and 
MD is resumed starting with the new configuration (see 



Fig-D- This ensures that the path generated from xq to 
Xp is fully reversible: a configuration in basin p reaching 
Xp would trigger the activation, bringing it to the other 
end of this path, at xo. Reversibility, in a weak sense, is 
ensured by the symmetric criterion for entering and leav- 
ing the saddle region as well as by keeping the path on 
a hypersurface of constant energy: for each transition, 
its inverse is also possible. In addition to reversibility, 
we have to ensure that the relative probabilities of these 
transitions are correctly weighted; this is discussed be- 
low. 

The conservation of energy requires that the length of 
the velocity vector as MD is restarted is equal to that at 
the beginning of the activation {\vp\ = |wo|); the direction 
should be chosen so as to point inside the new basin, but 
otherwise is arbitrary. After entering the new basin, MD 
is continued for a very small number of steps — 10 or 
so — to prevent the system from a quick recrossing back 
to the original basin. This is implemented by letting 
the system bounce back against the constant eigenvalue 
surface. After these few steps, the MD stage continues 
until the system crosses the basin boundary, and another 
activation is begun, etc. 

We note that the activation path does not always lead 
to a new basin. In our simulations, as is illustrated in 
Fig. n] it is not rare to see the trajectory form a circular 
path, coming back very close to the initial point xq, after 
a long excursion in the saddle region. It is also possible 
that the trajectory returns to the same basin but not at 
the starting point. This does not invalidate the algorithm 
but makes it less efficient. 



B. Calculating the event free energy 

Once we have a trajectory, it is necessary to compute 
the free energy difference between its beginning and its 
end. 

Consider the diagram in Fig. 12 This shows schemat- 
ically a few nearby activation paths between different 
basin boundaries (which are (TV — l)-dimensional hyper- 
surfaces). Points within area ^5*1 in the figure move to 
points within area dS2- If an ensemble with the micro- 
canonical distribution is considered, the density of flux 
of the trajectories through a hypersurface is the same 
for any hypersurface at all points having the same po- 
tential energy. Then the ratio of the rate of the direct 
transition (from dSi to dS'2) to the rate of the inverse 
transition (from ^5*2 to dSi) is equal to the ratio of the 
areas, dSi/dS2- If we want the microcanonical distribu- 
tion to be preserved, this ratio should equal 1. In general, 
however, it is not unity and we need to add an additional 
acceptance/rejection step, for instance, accepting a par- 
ticular activation transition with a Metropolis-like prob- 
ability Pace — min(c?52/dS'i, 1) for the transition from 
dSi to ^5*2 in the figure. 

The activation transition can be considered as a trans- 
formation between points on different basin boundaries 
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basin 1 




FIG. 2: Sketch of tube connecting two basins (see text). 



(for instance, point xq is transformed into point Xp, area 
dSi in Fig. 121 is transformed into area ^6*2). The ratio 
J = dS2/dSi corresponds therefore to the Jacobian of 
this transformation and it is all we need to ensure de- 
tailed balance. 



1. The Jacobian of the activation transformation 
boundary factor 



the 



Imagine again a tube formed by nearby activation 
trajectories (Fig. 121. If the trajectories forming the 
tube start within area dSi on the first basin boundary, 
then the cross-section at the beginning of the tube is 
dS'i = dSi cosai, where ai is the angle between the nor- 
mal to the initial basin boundary and the activation tra- 
jectory at its start. Likewise, the cross-section at the end 
of the tube is dS'2 — dS2 cos a2 , where a2 is the angle be- 
tween the normal to the final basin boundary and the 
activation trajectory at its end. Including all contribu- 
tions, the Jacobian can then be written as 



where 



Jb = 



cos ai 
cos a2 ' 
dS^ 
dS[' 



(3) 

(4) 
(5) 



We call Jb the boundary factor and Jxs the cross-section 
factor. 

We start by calculating the boundary part Jb- First, 
we note that Eq. ||2Jl for the activation trajectory is the 
discretized version of 



dx 
d7 



hix) + c{x)F{x). 



(6) 



This allows us to get an estimate of c, the factor preserv- 
ing the total energy during the activation. 
We write the change in potential energy as 

?^ = ?•VC/ = -?•F = -F||-cF^ (7) 
dr dr dr " 

where F|| = {F ■ h). Since we want to keep U constant, 
4^ = and 



El 

F2- 



(8) 



is a sum over N components of the force and thus 
scales as 0{N), the system size, since all modes are 
roughly equally excited. For its part, Fy represents just 
one component along the activated direction and does 
not grow with system size. These observations imply 
therefore that c scales as 0{1/N). 

Next, we will show that the eigenvector h is nearly 
tangent to the trajectory. Using Eq. 0, the angle (3 
between the activation trajectory and the eigenvector h 
is given by 



cos/3 



[h 



dx \ 



1 + cFi 



I ^ I 




Since, as discussed above, F^ ^ F^ 
cF\\,c'Fl 



(9) 

for big systems, 

<C 1, and cos/3 is nearly 1. Thus we can re- 
place the direction of the trajectory with the direction of 
h when calculating angles ai and a2- 

Note that since the basin boundary is by definition the 
constant-eigenvalue surface, the normal to it is parallel 
to VA, where A is the lowest eigenvalue. Then 



cos ai 2 



h ■ VA 



(10) 



where all quantities are evaluated at the beginning of the 
activation trajectory for ai and at its end for a2- Thus 
in order to calculate ai and a2, we need a way to find 
VA numerically. The most efficient method is as follows. 
By definition, at point x in the configuration space. 



H{x)h{x) X{x)h{x), 



(11) 



where H(x) is the Hessian operator at point x. Consid- 
ering X = {xi,X2, . . . , xat) as a set of parameters {xi}, 
we can apply the Hellmann-Feynman theorem and find, 
in vector form and with Einstein's summation conven- 
tion: 



VA hVHh 



dH^ ^ d^F 

oxi oxjOXk 



hjhk, (12) 



where Ci are unit vectors along the coordinate axes. 
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This expression is simply the second derivative of F 
along the direction of the eigenvector h with the negative 



sign, I.e. 



VA(f) 



lim 

(5^0 



2F{x) - F{x + (5 • /i) - F{x ~5-h) 
^2 • 



(13) 



It can be used directly for numerical evaluation. We sim- 
ply need to compute the force F at three nearby points 
along the direction of h for each boundary in order to 
obtain an accurate evaluation of the boundary factor Jb. 



Jacobian. However, straightforward evaluation of the di- 
vergences entering Eq. IjlTI) by calculating numerically 
the derivatives of h and F along N orthogonal directions 
for many points on the trajectory is computationally very 
costly and any usable method will require further careful 
analysis and making certain reasonable approximations, 
as discussed below. 

The logarithm of the cross-section factor in the Jaco- 
bian is an integral along the activation trajectory: 



In Jxs = 



j(x(r'))dr'. 



(18) 



2. Analysis of the cross-section Jacobian 



where 



The second factor in the total Jacobian is the cross- 
section factor Jxs. To evaluate it, we need to see how 
the cross-section of an infinitesimally narrow tube formed 
by activation trajectories changes between the two basin 
boundaries. Describing the evolution in the configuration 
space as a function of r by the equation 



dx - 

d7 = 



(14) 



then, as r is incremented by dr, point x transforms into 
X -\- f{x)dT. The Jacobian of that transformation is given 
by the determinant of the matrix Aij = Sij + ^-dr, 

which is l + T.i^dT + 0{dT^) ^l + d\Y fdr + OidT"^). 
For an infinitesimal volume 5V{x) around point af, the 
rate of change simply becomes 



—5Vix) = div/(f) • 5V{x). 
dr 



Which can be solved formally: 



5t/(T) = (5F(0) exp / div/(f(T'))dT' 



(15) 



(16) 



Going back to the continuous version of our evolution 
equation, Eq. ©, we note that \cF\ — ^ 1. 

Therefore, the speed along the activation trajectory is 
nearly constant (equal to one). Thus the infinitesimal 
volume 5V will only change its size and shape in the 
transverse directions, but will not shrink or expand in 
the longitudinal direction. Then the tube cross-section 
ratio between any two points on the trajectory is the 
same as the volume ratio between the same two points. 
The logarithm of the cross-section contribution to the 
Jacobian is then 



In X; 



In 



5V{t) 
5V{Q) ' JO 



div/(f(r'))dr' 



div K (x (r' ) ) + div c (x (r' ) ) {x{t' ) ) 



(17) 



dr'. 



Equation (|17|l . together with Eq. H13|l for the boundary 
factor, can be used in principle to calculate the activation 



j{x) = div ft. + div (cF) ^ div h + cdiv F + F ■ Vc. (19) 

Compare now the second and the third terms in Eq. 119|) 
to show that the third term can be neglected. In the 
second term, divF is the trace of the Hessian H taken 
with the negative sign and is therefore 0{N). Since c 
is 0{1/N), the second term in Eq. lO is 0(1). Now 
consider the third term. Using Eq. ((SJ, 



F-Vc = 



F^ 



F-VF^ F-VF, 



F2 



F2 



(20) 



If we use the coordinate system in which axes are parallel 
to the eigenvectors of the Hessian at point x (in partic- 
ular, the zeroth axis is parallel to h), then dFi/dxj = 
—XiSij, where Ai is the ith eigenvalue of the Hessian. 
Then the first term in Eq. if^ is 



F\\ F-VF^ 



F2 



= 2c 



(21) 



which is 0{\/N) (given that aU A's are 0(1), c is 0{1/N) 

and Y^^^a^ — ^^"^ thus negligible compared 
to the second term of Eq. (|19|l . In the second term of 
Eq. 



F ■ Vi^ii 



F- V 



F,h, 



F2 

F-j:Z='o{im/dx,)h,e,) 

J J. 

F2 

y-W-i din p.p. 



p2 



p2 



(22) 



In the last expression, the first term is clearly 0{1/N) 
and thus negligible; the second term is also negligible 
(this will be so even under a completely unrealistic as- 
sumption that all of dhi/dxj are 0(1), since Fi are of 
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random signs). Thus the third term in Eq. (|19|l can al- 
ways be neglected for big enough N and we end up with 

j divK + cdivF. (23) 

In Appendix A, we will discuss the physical meaning 
of the second term in Eq. (|23|) . using the harmonic ap- 
proximation. 

C. Implementation of the POP-ART algorithm 

The actual implementation of POP-ART, as used to 
obtain the results presented in the next section, incorpo- 
rates the following steps: 

1. We start with MD at finite temperature and first 
equilibrate by rescaling the velocities. We use a 1 
fs step and compute the lowest eigenvalue every 10 
steps. After we have crossed the basin boundary 
defined by the threshold, we retrace our MD path 
and identify the crossing time with an accuracy of 
1 fs. 

2. We then apply Eq. ^ and generate the event from 
one basin to another, stopping at the threshold and 
saving configurations along the way. We take At = 
0.01 A. 

3. From the first and final configurations of the acti- 
vation path, we compute the boundary factor, Jb, 
using Eqs. JTHl and IpH) . 

4. We then evaluate the cross-section Jacobian Jxs by 
integrating j, as defined by Eq. (|23|) . over the path- 
way. We now have the full free energy difference 
between the entry and exit points. 

Straightforward evaluation of the divergences en- 
tering Eq. H23|l by calculating numerically the 
derivatives of h and F along N orthogonal direc- 
tions for many points on the trajectory is compu- 
tationally very demanding. It is possible, how- 
ever, to lower this cost significantly while keep- 
ing a reasonable accuracy. First, we use a 15- 
iteration Lanczos scheme which allows us to ob- 
tain the eigenvector h within 0{N). The diver- 
gence of the eigenvector can also be obtained with 
0{N), admittedly with a much larger prefactor, 
provided the potential is sufficiently short-ranged. 
For two atomic coordinates i and j belonging to 
atoms which are well outside of each other's inter- 
action range, dhi/dxj « 0. We can exploit this 

property, to obtain divh = ^^^q^ dhi/dxi with 
less than 0{N) force evaluations. For example, 
two terms in this summation can be obtained with 
one eigenvector computation: dhi/dxi+dhj/dxj w 
T,m={i,j}ihmix + Aci + Ae,) - hm{x))/A. This 
trick can easily be extended as long as the added 
coordinates are sufficiently far apart. For the two 



systems studied here, the unit cell is divided into 25 
groups of 40 non-interacting atoms each. The to- 
tal cost of evaluating div h adds up to 75 Lanczos 
recursions with only 5 iterations each. 

5. The previous two steps provide Jb and Jxs and thus 
the free energy difference between the first and the 
last state on the activation trajectory, which is then 
used in a Metropolis accept/reject move. If the 
event is rejected, the component of the velocity nor- 
mal to the basin boundary is reversed, and MD is 
continued. If the event is accepted, we continue 
MD in the new basin, using the initial velocity (re- 
versing the component of the velocity normal to 
the basin boundary, if necessary). As mentioned 
above, we run 10 steps to bring the configuration 
away from the border. 

6. Once the threshold is reached again, repeat steps 
2-6. 



IV. SIMULATION RESULTS 

To verify the thermodynamical correctness of the 
POP-ART method and to investigate its computational 
efficiency, we have studied the diffusion by interstitials 
and vacancies in a silicon crystal, described by the 
Stillinger- Weber potential 13j in systems of respectively 
1001 and 999 atoms, and a cubic simulation cell of 
27.136 A3. 



A. Interstitial diffusion 

First, we look at the interstitial in a 1001-atom cell 
of Si. In principle, many interstitial sites are possible in 
this system, but only three of them are significantly pop- 
ulated: the lowest-energy configuration and two config- 
urations with almost the same ener gy ( within 0.01 eV), 
about 0.75 eV above the first one |2y|- We are inter- 
ested in computing the probability of being in one of the 
higher-energy states. 

If the population of the higher-energy states were 
determined by the energy difference alone, this would 
amount to a population of the higher-energy states of 
only 0.07% at 1200 K. However, there are degeneracies 
and the potential wells of the higher-energy states are 
much flatter than that for the low-energy state, lead- 
ing to a noticeable entropy difference. Because of this 
difficulty, we extract the thermodynamical equilibrium 
between these two states with MD. This forces us to per- 
form the tests at a relatively high temperature. Here, we 
report results for 1000 K and 1200 K. 

These conditions are not ideal for POP-ART since at 
such high temperatures the jumps between the minima 
are rather frequent and straightforward MD is quite effi- 
cient. But given the significance of both energetic and en- 
tropic contributions, as well as some anharmonic effects 



7 



present at such high temperatures, it is a very good test 
for the accuracy (rather than efSciency) of POP-ART. 

At 1000 K and 1200 K, the system can spend a non- 
neghgible amount of time outside the basins, i.e., in the 
saddle region (which would not be the case in systems 
more appropriate for POP-ART), leading to a different 
value of the ratio of time spent in the upper vs. the lower- 
energy states. We therefore need to distinguish carefully 
between the probability of being in the attraction region 
of a given minimum and the time spent in a particular 
basin (understood as defined in this paper) as a fraction 
of the total time spent within all basins. It is the latter 
quantity that we use to compare to the POP-ART result. 

At 1000 and 1200 K, the MD result for this quantity 
is, respectively, 1.6±0.1% and 3.6±0.1%, determined as 
an average over 25 runs, each of which lasted 10 ns and 
an assignment to a basin, with threshold X — —2 and 
— 1 eV/A^, done every 100 fs. 

For POP-ART the fraction is obtained as an average 
over 5 runs, of 10000 iterations of the POP-ART algo- 
rithm each. The assignment to a basin is done every 100 
fs, excluding the 10 steps for moving away from the bor- 
der after each generated event. At 1000 and 1200 K, we 
obtain the respective ratios of 1.3 ± 0.3% and 3.5 ± 0.3%. 
Clearly POP-ART samples accurately the thermodynam- 
ical weight of the high and low-energy local mimina. 
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FIG. 3: Diffusion of a vacancy in Stillinger- Weber silicon 
(counted as the number of jumps per million force evalua- 
tions) as a function of temperature for molecular dynamics 
and POP-ART. The dashed line is an Arrhenius fit to the 
MD results with an energy barrier of 0.46 eV (close to the 
value of 0.43 eV reported in Ref. 12^1. and the solid line is a 
1/T fit to the POP-ART results. 



to the fact that before POP-ART attempts to find its way 
to a new basin, it needs to reach a given curvature thresh- 
old. At a very low temperature, even this will be a very 
rare event. 



B. Vacancy diffusion 

Having established the accuracy of POP-ART, wc now 
characterize its efficiency. For this, we consider a vacancy 
in a 999-atom cell of Si. Vacancy diffusion is associated 
with a single activation barrier of 0.43 eV [l?!]. Assessing 
the efficiency of POP-ART relative to MD can be cleanly 
done in this system, since the speed of phase space ex- 
ploration is simply given by the vacancy's diffusion coef- 
ficient, which follows the Arrhenius law. The comparison 
with MD is done on the basis of the number of calls to 
the force routine, since that takes ~99% of the computer 
time. 

Figure 131 shows an Arrhenius plot of the diffusion rate 
per million force operations obtained by MD and by 
POP-ART, as a function of temperature. The diffusion 
rate is defined as the vacancy hopping rate, which is lin- 
early related to the square displacement per unit time. 
The value of A used in the POP-ART simulation is se- 
lected to optimize the diffusion rate. The threshold val- 
ues used in the simulation are A = — 2 eV/A^ for T > 
750 K, A = -1 eV/A2 for 600 K < T < 750 K and 
A = eV/A2 for T < 600 K. Clearly, the diffusion rate 
per force evaluation obtained with POP-ART does not 
show activated behavior and provides a significant boost 
at temperatures below 700 K, reaching a factor of more 
than 4 orders of magnitude at room temperature. Inter- 
estingly, the diffusion rate per force evaluation is not con- 
stant with POP-ART but can be fitted by a 1/T curve. 
The slowing down with decreasing temperature is related 



V. DISCUSSION AND CONCLUSION 

Efficient sampling of slow systems is one of the main 
challenges in computational physics today. For disor- 
dered systems such as glasses, for example, standard 
thermodynamical methods fail because of the very small 
phase space occupied by the relevant configurations. Ac- 
tivated methods, such as ART, overcome these difficulties 
by generating physically-possible trajectories through the 
conformation space but they do not offer a proper ther- 
modynamical sampling. 

The properly-obeying-probability activation- 
relaxation technique (POP-ART) lifts these limitations 
by generating activated trajectories with proper thermo- 
dynamical weighting. Mixing molecular dynamics with 
activation over barriers, this algorithm respects detailed 
balance and samples in a well-defined thermodynamical 
ensemble. 

To verify the correctness of POP-ART, we sampled the 
states visited by an interstitial in c-Si. In this system, 
the higher-energy states of the interstitial are energet- 
ically suppressed but entropically favored. Comparing 
with MD, we found that POP-ART samples the high- 
energy states with the proper probability, demonstrating 
its correctness. In order to assess the efficiency of this 
method, we looked at the vacancy diffusion also in c-Si. 
In this case, POP-ART is found to outperform MD below 
700 K, and it is about four orders of magnitude faster at 
room temperature. 
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One of the advantages of POP-ART is that all the 
information it needs is local. This makes it possible to 
apply a number of approximations to increase further its 
efficiency without sacrificing the sampling. POP-ART 
can also be extended to reproduce the correct activated 
dynamics; this is currently examined and will be reported 
in a further publication. 
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Appendix A: The physical meaning of the 
divergence of the force 

Using the harmonic approximation, we can reveal bet- 
ter the physical meaning of the second term in Eq. (|23|l . 
First of all, using the expression for c, Eq. we get 



cdivi^ 



F2 



divF. 



(24) 



At relatively low temperatures, the system is well de- 
scribed by the harmonic approximation. The potential 
around the minimum can then be rewritten as 



1 



(25) 



1=0 



where Xi represent the normal modes, and the force be- 
comes 



divF w - ^ fci 



(26) 



Using the same approximation for and neglecting an- 
harmonicity, we get: 



2 2 
i ■ 



(27) 



i=0 



Since the spectrum should be dense, xf can be replaced 
with their thermal averages {xD — kBT/ki, giving 



N-l 



and 



cdivi^ 



i=0 



keT 



(28) 



(29) 



After integrating over the whole activation trajectory, we 
obtain 



zdYvFdr' = - 



knT' 



where 



F\\dT'. 



(30) 



(31) 



Then the contribution of the c div F term to the Jacobian 
J is 



exp(|cdivFdr')=exp(-0), 



(32) 



a Boltzmann factor. Ai?|| is essentially the energy change 
that would have occurred along the activation trajec- 
tory, if it were parallel to h everywhere and the energy- 
correcting cF term in Eq. ^ was not there. 

The transition probability between two minima should 
contain both energetic and entropic contributions. Given 
the above result, it is tempting to associate the div h term 
with the entropic and the c div F with the energetic con- 
tribution; however, an example considered in Appendix B 
shows that the reality is more complex: in fact, Ai?||, de- 
fined as above, is temperature-dependent and the cdivi*" 
term therefore contains both energetic and the entropic 
contributions. 



Appendix B: the Jacobian in a model potential 

To get some insight into the physical meaning of the ac- 
tivation Jacobian and its different components, consider 
the following model, defined by the potential: 



N-l 



U = Uo{xo) + ^ ^ h{xo) 



(33) 



Here Uq{xq) is a function with two minima and a maxi- 
mum between them, so that coordinate xq describes the 
activated mode (i.e., is the "reaction coordinate"), and 
the other A^ — 1 degrees of freedom serve as the "heat 
bath". "Force constants" ki{xo) are assumed to remain 
positive for all xq of interest (e.g., between the minima). 

Our model does not represent the most general situ- 
ation. In particular, the eigenmode-foUowing transition 
path (such as, e.g., ART would find) is a straight line (co- 
inciding with the 0th axis); also, on that line, an eigenvec- 
tor for a particular mode has the same direction (parallel 
to a coordinate axis) everywhere. However, the model is 
interesting enough: as the frequencies of the bath modes 
(determined by ki(xo)) can change along the activation 
path, there are both energetic and entropic contributions 
to the probability of being in a particular place along the 
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reaction coordinate. Indeed, the probability density of 
having the zeroth coordinate equal to xg is 

p{xq) cx J c?xi . . . dxAT-i exp[— C///cb2^] 
UoixoY 



exp 



X / dxi . . . dxN-i exp 



N-l 

E 



ki{xo)xf 



2kBT 



exp 



Uo{xo) 



knT 



or 



p{xq) oc exp[-T{xa)/kBT], 
where the free energy 

T{xq) = Uo{xo) - TS{xo) 
and the entropy 



(34) 

(35) 
(36) 

(37) 



We will assume in what follows that ki{xQ) are linear 
functions, i.e.. 



(38) 



The matrix elements of the Hessian for the potential 
given by Eq. are 



i^oo = C/o(^o), 
Ha = ki{xo), i ^ 0, 



(39) 



and the remaining elements are zero. The force compo- 
nents are 



AT-l 
i=l 

-fcj(a;o)a;i, i ^ 



(1)^2 



F, = -h{xo)x„ 1^0. (40) 
For a Hessian with only Ha and iJoi = Hiq non-zero. 



h = C[l, 



22 



(41) 



where C is the normalization constant, and the eigen- 
value A is the lowest solution of the following equation: 



Af-l 



i=l 



A — Hi, 



(42) 



In a real physical system, the entropy change along 
an activation trajectory will always remain finite and of 



the same order of magnitude as the thermal energy, as 
the system size increases. In our model, this will mean 

that most fcj"'^'' are small enough. Note that this is essen- 
tially the same as the assumption of most modes being 
nearly harmonic that we have used when approximating 
the c div F term in the Jacobian. Specifically, 



^1 fc(i) 



TAS - ksT J2 ^^2;o - fc^T, 



or 



^1 fc(i) 



E ^A.o - 1. 



k 

For simplicity, we will assume in addition that 



E 



Hot 



1=1 

7V-1 



^^00 — Ffi, 



E 



1 {Hqo - Hi; 



^ Hqo, 



< 1. 



(43) 



(44) 



(45) 



(46) 



This will be the case, in particular, at low enough T, 
when the magnitudes of most small. Under these 

conditions, 



A « i/oo - U^lixo), 
C « 1. 



(47) 
(48) 



The constant-eigenvalue surfaces are then orthogonal to 
the 0th axis and the cosine of the angle between the nor- 
mal to a constant-eigenvalue surface and h is nearly 1, 
so the boundary contribution to the activation Jacobian 
can be neglected. Consider now the cross-section factor. 
We need to calculate j, as given by Eq. if^ . First of all. 



k['^x. 



Ui;ixo)~ki{xo)' 



and so 



E 77777 



.(1) 



^ U(l{xo) - ki{xo) ' 
Further, using Eq. 

F|| ^ {F-h) = -Ul,{xo)--Y.ki'^x: 



(49) 



(50) 



^ ki{xo)k^^^xf 
Ul{{xo)~k,{xo) 

We can now use Eq. H29|) to calculate cdivF; it is, how- 
ever, instructive to repeat the considerations, to see how 
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exactly the appropriate approximations are made in this 
particular case. We obtain 



N-1 



(where an approximation similar to Eq. H26|l is obtained 
by neglecting the term associated with the anharmonic 
activated mode) and 



N-l 



N-1 



(53) 



In the last expression, bearing in mind the condition (|44|l . 
the second term in the parentheses is of the same order 
of magnitude as the first one, and then both of them 
are negligible compared to the second sum (this is an 
approximation similar to the one used in Eq. (|27|l l. We 
therefore obtain 



cdivF 



C^o(-o) + iE^ 



^o'(^o) + hjxo) 



1 Uaixo) ~ h{xQ) 



(54) 



As explained in Appendix A, we can replace with their 
thermal averages and then 



[Uoixo) 

l^V^'(xo) + fe.(xo) kl 



(55) 
'keT. 



Then finally 



j = -Uo{xQ)/kBT - i 



2 ki 



(56) 



and the Jacobian J is 



J = exp{Jj{xo)dxo)^exp[-{AUo-TAS)/kBT] 
= eyii>[-AT/kBT], (57) 

where AUq is the change in Ua, AS is the change in en- 
tropy S [as given by Eq. ifTfjl ] between the two basin 
boundaries, and AT = AUq — TAS — exactly as ex- 
pected. 
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